!	The parameters estimated are >>>out-of-sample parameters<<< 
!	the difference of this code from brand_upd2 is that subroutine used in this program
!	the out-of-sample parameters are estimated by year, this code is verfied
!   2/20/13 ::  New issues:  1) Use a minimum of 15 years of ``in-sample'' data (to handle the over-fitting problem).
!                            2) Line up the dates properly.
!
!   2/11/19 ::  New issues:  1) I now have the permnos in the input files.
!                            1) I also have an SIC indicator for financials in the last column.
!   3.31.22 ::  I am starting a new project to look at the nature of overfitting.  First step is to run the optimal mode
!               for the gamma 2 investor using gamma* = 3 to see if things change with 36 more months of data.
program main
USE SHARE
USE RNUND_INT
USE UMIAH_INT
! USE UVMID_INT
IMPLICIT NONE
INTEGER :: dumfnm,ik,ip,lp,i,j,errcod, STATUS, iostat,T,K,P,rednfm,NBS,ib,js,ic
INTEGER :: IPARAM(7),dumper, jm,dumnfm,NP,jp,dnsm,ljm,ja,NYP,INTERK,NNTRK,permnd,finind
DOUBLE PRECISION, DIMENSION(5000,900,10) :: ochars
DOUBLE PRECISION, DIMENSION(5000,900) :: oweig,orit
DOUBLE PRECISION, DIMENSION(900,10) :: chrsum,sschr
DOUBLE PRECISION, DIMENSION(900) :: sumsz
DOUBLE PRECISION :: rdum,momdum,btmdum,medum,sew,sow,svw,ssow,ssew,ssvw,sw,dumome,lmdum
DOUBLE PRECISION:: sumpar(10),sspar(10),sumc(10),ssumc(10)
DOUBLE PRECISION:: U,betdum,sbdum,FUN,FOU,sbsw,FX,rsddum,lszdum,BMdum,LYSDUM,afydum
DOUBLE PRECISION::para(900,10),EU(900),ttu, optu, lastu
DOUBLE PRECISION:: opw(5000,900), opret(900), mpret(900), vpret(900)
DOUBLE PRECISION:: H(10,10),G(10),THETA(10)
DOUBLE PRECISION:: dbsumw,snwtp,avgnw
INTEGER, DIMENSION(5000) :: IR
EXTERNAL FCN,GRAD,HESS
OPEN(unit= 14, file='/Users/chrislamoureux/BSV/22/data/BSVUSE22.dat', STATUS ='OLD',iostat= errcod)
OPEN(unit= 17, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M7/g9/WTS.dat', STATUS='replace', iostat= errcod)
OPEN(unit= 21, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M7/g9/theta.dat', STATUS='replace', iostat= errcod)
OPEN(unit= 22, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M7/g9/rets.dat', STATUS='replace', iostat = errcod)
OPEN(unit= 23, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M7/g9/nstks.csv', STATUS='replace', iostat = errcod)
OPEN(unit= 33, file='/Users/chrislamoureux/BSV/22/data/output/bootstrap/M7/g9/ISrets.csv', STATUS='replace', iostat = errcod)
CALL ERSET(3,0,-1)
zero=0.0d0
one=1.0d0
two=2.0d0
NP=5         ! number of parameters (i.e., length of theta vector)
gamma=9.0d0  ! gamma* (i.e., the concavity of the loss function used to generate the feasible portfolio rule).
T=744   ! number of months in total.  This does not include the first data collecting period.
        !  i.e., here the sample runs from month 61 through 804 (through 2021).
P=180	! The next month (P+1) is the starting month for out-of-sample portfolios.
NYP=(T-P)/12
!
DO i=1,T
   ktr(i)=0
   DO j=1,NP
      chrsum(i,j)=zero
      sschr(i,j)=zero
   END DO
   sumsz(i)=zero
END DO
!
ljm=0
DO
   READ (14,2166,iostat=status) jm,rdum,momdum,BMdum,lszdum,betdum,lysdum,afydum,rsddum,finind,dumome
            2166 format (i3,1x,8(f12.6,1x),i1,1x,f17.4)
   IF(status /= 0)  EXIT
   IF (jm /= ljm) THEN
      ktr(jm)=1
   ELSE
      ktr(jm)=ktr(jm)+1
   END IF
   orit(ktr(jm),jm) = rdum    !  read one stock sequentially
   oweig(ktr(jm),jm)= dumome
   sumsz(jm)=sumsz(jm)+dumome
   ochars(ktr(jm),jm,1) = momdum              ! 1
   ochars(ktr(jm),jm,2) = BMdum               ! 2
   ochars(ktr(jm),jm,3)  = lszdum             ! 3
!    ochars(ktr(jm),jm,2)  = betdum             ! 4
!*Off*****    ochars(ktr(jm),jm,1)  = lysdum
   ochars(ktr(jm),jm,4)  = afydum             ! 5
   ochars(ktr(jm),jm,5)  = rsddum             ! 6
   DO j=1,NP
      chrsum(jm,j)=chrsum(jm,j)+ochars(ktr(jm),jm,j)
   END DO
   ljm=jm
END DO
DO i=1,T
   PRINT 1126, i,ktr(i); 1126 format (2(2x,i4))
   WRITE (23,1126) i,ktr(i)
END DO
! Standardize and normalize the parameters and define the market weights.
DO i=1,T
   DO j=1,NP
      chrsum(i,j)=chrsum(i,j)/ktr(i)
   END DO
   do j=1,ktr(i)
      oweig(j,i)=oweig(j,i)/sumsz(i)
      WRITE (17,1766) i,j,oweig(j,i); 1766 format (i3,2x,i4,f12.9)
   end do
END DO
DO i=1,T
   DO j=1,ktr(i)
      DO ja=1,NP
         sschr(i,ja)=sschr(i,ja)+(ochars(j,i,ja)-chrsum(i,ja))*(ochars(j,i,ja)-chrsum(i,ja))
      END DO
   END DO
END DO
DO i=1,T
   DO j=1,NP
      sschr(i,j)=sqrt(sschr(i,j)/(ktr(i)-one))
   END DO
END DO
DO i=1,T
   DO j=1,ktr(i)
      DO ja=1,NP
         ochars(j,i,ja)=(ochars(j,i,ja)-chrsum(i,ja))/sschr(i,ja)
      END DO
   END DO
END DO
      
!
!   out of sample series
!   IF (T /= 100) STOP 1091
NBS=10000
INTERK=0
NNTRK=0
DO ib=1,NBS
   INTERK=INTERK+1
   IF (INTERK == 100) THEN
      NNTRK=NNTRK+1
      INTERK=0
      PRINT 789, NNTRK; 789 format (' First ',i3' Bootstrap samples finished.')
   END IF
! Construct the pseudo sample for this bootstrap:
   
   DO jm=1,T
      DO ic=1,NP
         sumc(ic)=zero
         ssumc(ic)=zero
      END DO
      sbsw=zero
      CALL RNUND(ktr(jm),ir)
      DO js=1,ktr(jm)
! For debugging purposes use the sample:
!
!          ir(js)=js
!
         rit(js,jm)=orit(ir(js),jm)
         weig(js,jm)=oweig(ir(js),jm)
         sbsw=sbsw+weig(js,jm)
         DO ic=1,NP
            chars(js,jm,ic)=ochars(ir(js),jm,ic)
            sumc(ic)=sumc(ic)+chars(js,jm,ic)
         END DO
      END DO
      DO ic=1,NP
         sumc(ic)=sumc(ic)/ktr(jm)
      END DO
      DO js=1,ktr(jm)
         weig(js,jm)=weig(js,jm)/sbsw
         DO ic=1,NP
            ssumc(ic)=ssumc(ic)+((chars(js,jm,ic)-sumc(ic))*(chars(js,jm,ic)-sumc(ic)))
         END DO
      END DO
      DO ic=1,NP
         ssumc(ic)=sqrt(ssumc(ic)/ktr(jm))
      END DO
      DO js=1,ktr(jm)
         DO ic=1,NP
            chars(js,jm,ic)=(chars(js,jm,ic)-sumc(ic))/ssumc(ic)
         END DO
      END DO
   END DO
   DO jp=1,NP
     theta(jp)=zero
     sumpar(jp)=zero
     sspar(jp)=zero
   END DO
   DO k=1,NYP           ! We update theta at the beginning of each year using all historical data.
      NMOS=P+(k-1)*12   ! NMOS is the in-sample period.
!       PRINT 4266, NMOS,P; 4266 format (' in main before umiah nmos: ',i4,' p ',i4)
      IPARAM(1)=0
      CALL UMIAH (FCN,GRAD,HESS,theta,N=NP,IPARAM=IPARAM,FVALUE=FOU)   ! IMSL Math V2 p. 1213.
      WRITE (33,3366) NMOS; 3366 format (i4)
      DO i=1,NMOS
         WRITE (33,3367) i,ISOPRT(i); 3367 format (i4,1x,f14.10)
      END DO
!         CALL UVMID(FCN,GRAD,-12.d0,22.d0,theta,FX=FX)
!       PRINT 1426, FOU,(theta(jp),jp=1,NP),IPARAM(3)
            1426 format (' UMIAH finished Max U: ',f12.7,' theta: ',5f12.7,' no iterations: ',i4)
           FOU=FX
           EU(k)=FOU
   	   WRITE(21,114) k, (theta(jp),jp=1,NP); 114 FORMAT(1X, I8, 10(1X, F18.8))
         DO i = 1, 12	! assign parameters to the next year
            DO jp=1,NP
               para(NMOS+i,jp) = theta(jp)
               sumpar(jp)=sumpar(jp)+ theta(jp)
            END DO
         END DO
   END DO   ! In-sample utility optimization for first NMOS months complete.
!
!    Now set up the next 12 months -- the oos period corresponding to the k'th optimization:
!
   U=zero
   avgnw=zero
   DO i=P+1,T
      snwtp=zero
      sw=zero
      DO j=1,ktr(i) ! optimal weights and returns of stocks
         sc=zero
         DO jp=1,NP
            sc=sc+para(i,jp)*chars(j,i,jp)
         END DO 
         opw(j,i)=weig(j,i)+sc/ktr(i)
         IF (opw(j,i) < zero) THEN
            snwtp=snwtp+opw(j,i)
         END IF
         sw=sw+opw(j,i)
!        WRITE(20,121) MDATE(i),ktr(i),weig(j,i),opw(j,i); 121 FORMAT(1X,I8,1X, I5,2(2X,F16.13))
      END DO
      avgnw=avgnw+snwtp
!  PRINT 5632, i,MDATE(i),sw; 5632 format (' Month ',i4,' date: ',i8,' sum wts: ',f12.10)
      opret(i)=zero; mpret(i)=zero; vpret(i)=zero
      dbsumw=zero
      DO j =1,ktr(i) ! Portfolio returns in each month
         opret(i)=opret(i) + opw(j,i)*rit(j,i)
         mpret(i)=mpret(i) + rit(j,i)
         vpret(i)=vpret(i) + weig(j,i)*rit(j,i)
         dbsumw=dbsumw+opw(j,i)
      END DO
      PRINT 8761, i,dbsumw; 8761 format (2x,i4,' debugging sum of opt wts: ',f12.9)
      U=U +((one+opret(i))**(one-gamma))/(one-gamma)
   END DO
   U=U/(T-P)
   avgnw=avgnw/(T-P)
   PRINT 5439, U; 5439 format (' Mean oos utility: ',f12.9)
   PRINT 5446, avgnw; 5446 format (' Avg neg wts: ',g16.7)
   WRITE (17,1776) avgnw; 1776 format (f16.8)
   sew=zero
   svw=zero
   sow=zero
   ssew=zero
   ssvw=zero
   ssow=zero
!    DO i=1,T
!       PRINT 7651, i,MDATE(i),ktr(i); 7651 format (1x,i4,2x,i8,2x,i7)
!    END DO
   DO i = P+1,T
      mpret(i)=mpret(i)/ktr(i)
      sew=sew+mpret(i)
      svw=svw+vpret(i)
      sow=sow+opret(i)
      WRITE(22,117) i,opret(i),mpret(i),vpret(i); 117 FORMAT(1X,I8,4(2X,F12.8))
   END DO
   sew=sew/(T-P)
   svw=svw/(T-P)
   sow=sow/(T-P)
   DO i=1,NP
      SUMPAR(i)=SUMPAR(i)/(T-P)
   END DO
   DO i=P+1,T
      ssew=ssew+(mpret(i)-sew)*(mpret(i)-sew)
      ssvw=ssvw+(vpret(i)-svw)*(vpret(i)-svw)
      ssow=ssow+(opret(i)-sow)*(opret(i)-sow)
      DO jp=1,NP
         sspar(jp)=sspar(jp)+(para(i,jp)-sumpar(jp))*(para(i,jp)-sumpar(jp))
      END DO
   END DO
   DO jp=1,NP
      sspar(jp)=sqrt(sspar(jp)/(T-P))
   END DO
   ssew=sqrt(ssew/(T-P))
   ssvw=sqrt(ssvw/(T-P))
   ssow=sqrt(ssow/(T-P))
!   PRINT 7861, gamma; 7861 format (' Gamma ',f8.3)
   PRINT 6574, ib; 6574 format (' Bootstrap sample: ',i6)
   PRINT 7865; 7865 format (28x,' Mean ',8x,'Standard Deviation')
   PRINT 7866, sew,ssew; 7866 format(' Equally weighted Index',2(4x,f12.9) )
   PRINT 7867, svw,ssvw; 7867 format(' Value weighted Index  ',2(4x,f12.9) )
   PRINT 7868, sow,ssow; 7868 format(' BSV Optimal Portfolio ',2(4x,f12.9) )
   DO jp=1,NP
      PRINT 9436, jp,sumpar(jp),sspar(jp); 9436 format (1x,i2,2(2x,f13.8))
   END DO
! 1 bootstrap sample finished.
END DO
STOP; END
